Time Series Project

Spring 2024

Author

Nicholas Sager, Anishka Peter

Published

April 8, 2024

Data

Introduce dataset, response variable, and scenario. ACFs and Spectral Density plots.

library(dplyr)
library(lubridate)
bike_data <- read_csv("https://raw.githubusercontent.com/NickSager/MSDS-6373-Time-Series/master/Project/bike_data.csv")
Rows: 17379 Columns: 14
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): Date
dbl (13): Season, Hour, Holiday, Day of the Week, Working Day, Weather Type,...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
plotts.sample.wge(bike_data$`Total Users`)$xbar

[1] 189.4631
# make the bike data daily instead of hourly 
daily_bike_data <- bike_data %>%
  group_by(Date = as.Date(Date, format = "%m/%d/%Y")) %>%
  summarise(
    Season = mean(Season),
    Hour = mean(Hour),
    Holiday = mean(Holiday),
    Day_of_the_Week = mean(`Day of the Week`),
    Working_Day = mean(`Working Day`),
    Weather_Type = mean(`Weather Type`),
    Temperature = mean(`Temperature F`),
    Temperature_Feels = mean(`Temperature Feels F`),
    Humidity = mean(Humidity),
    Wind_Speed = mean(`Wind Speed`),
    Casual_Users = sum(`Casual Users`),
    Registered_Users = sum(`Registered Users`),
    Total_Users = sum(`Total Users`), 
  )

plotts.sample.wge(daily_bike_data$Total_Users)$xbar

[1] 4504.349

Models

Univariate: At least 2 candidate ARMA / ARIMA models

Seasonal Model

# create model 
plotts.sample.wge(daily_bike_data$Total_Users)$xbar # looking at the ACF's there seems to be a seasonality of 2

[1] 4504.349
factor.wge(phi = c(0,1)) 
  
  
Coefficients of AR polynomial:  
0.0000 1.0000 

                           AR Factor Table 
Factor                 Roots                Abs Recip    System Freq 
1+1.0000B             -1.0000               1.0000       0.5000
1-1.0000B              1.0000               1.0000       0.0000
  
  
est.arma.wge(daily_bike_data$Total_Users, p = 15) # looking at the overfit tables there's evidence of seasonality of 2 
  
  
Coefficients of AR polynomial:  
0.4739 0.0471 0.0189 0.0204 0.0573 0.1365 0.0094 0.0292 -0.0067 0.0560 -0.0099 0.0108 -0.0145 0.1114 0.0367 

                           AR Factor Table 
Factor                 Roots                Abs Recip    System Freq 
1-0.9950B              1.0050               0.9950       0.0000
1-1.1706B+0.7891B^2    0.7417+-0.8468i      0.8883       0.1355
1-0.5431B+0.7505B^2    0.3618+-1.0961i      0.8663       0.1993
1-1.5898B+0.7490B^2    1.0613+-0.4570i      0.8654       0.0647
1+0.8530B             -1.1723               0.8530       0.5000
1+0.3080B+0.7063B^2   -0.2180+-1.1698i      0.8404       0.2793
1+0.9479B+0.6822B^2   -0.6947+-0.9916i      0.8260       0.3473
1+1.4043B+0.6498B^2   -1.0805+-0.6094i      0.8061       0.4183
1+0.3114B             -3.2113               0.3114       0.5000
  
  
$phi
 [1]  0.473930982  0.047079867  0.018920349  0.020428876  0.057306483
 [6]  0.136485202  0.009418334  0.029186822 -0.006667372  0.055969333
[11] -0.009932710  0.010806592 -0.014549844  0.111431669  0.036710926

$theta
[1] 0

$res
  [1] -4.355839e+02 -4.798747e+02  1.600982e+02  1.312135e+02  6.391881e+01
  [6]  7.281801e+01  1.311062e+01 -4.872231e+02 -4.320116e+02  1.297026e+02
 [11] -1.614520e+02 -2.367857e+02  6.806097e+01  5.340858e+01 -1.089157e+02
 [16] -8.544696e+01 -3.025770e+02 -5.464649e+02  5.478171e+02  3.675432e+02
 [21] -1.527562e+02 -4.744995e+02 -1.297731e+02  2.670819e+02  5.051950e+02
 [26] -1.251830e+03 -6.196546e+02  2.834141e+02 -1.584651e+02 -2.024688e+02
 [31]  2.270164e+02  1.502517e+02  2.399052e+02  4.552043e+01  1.704674e+02
 [36] -4.924939e+02  4.168284e+02  1.822961e+02 -1.841704e+02  1.027864e+02
 [41]  5.863149e+01  2.876028e+02 -2.012276e+02  6.008407e+01  2.784402e+02
 [46]  6.306895e+01  3.438165e+02  5.445570e+02  8.139179e+02 -6.768341e+02
 [51]  1.264121e+01 -7.918039e+02 -1.409036e+02  1.327638e+02 -2.117107e+02
 [56] -3.834864e+02  2.987822e+02  6.011032e+02 -7.227381e+02  1.498051e+02
 [61]  2.308145e+02 -3.437019e+02 -6.842273e+01  6.444707e+01 -1.369469e+03
 [66]  6.117052e+02  3.115837e+02 -1.447301e+02 -1.306133e+03  7.212810e+02
 [71]  3.709251e+02  3.142353e+02 -1.141731e+02  7.815020e+01  3.448711e+02
 [76]  6.166459e+02  8.396201e+02  3.438872e+02  1.438669e+01 -2.704115e+02
 [81]  4.667250e+02 -5.165525e+02 -3.834151e+02  3.092290e+01  1.984874e+02
 [86] -7.812180e+02 -1.034169e+02  2.758276e+02 -7.858865e+02 -3.064390e+02
 [91]  9.060128e+01 -3.457090e+01  9.459615e+02  4.050910e+02 -8.895343e+02
 [96]  7.658151e+02  6.399149e+02 -1.318911e+03  2.990799e+02  4.741825e+02
[101]  7.702361e+02 -9.651203e+02 -1.605384e+02  1.099190e+03  2.935968e+02
[106] -2.119036e+03  1.779012e+03  3.626393e+02  2.078548e+02  9.021908e+02
[111]  7.525012e+02 -1.513207e+03  1.588468e+03  6.968067e+02  2.583676e+02
[116]  8.053814e+02  4.799753e+01  5.919175e+02  6.720277e+02  1.434884e+03
[121] -1.224971e+03  7.628568e+02  3.339877e+02 -1.625062e+03  8.476300e+02
[126]  5.166007e+02  5.407134e+02 -1.652893e+02  1.384048e+02  6.022095e+02
[131] -2.965421e+02  5.932004e+02 -5.550719e+02 -9.017327e+02  6.682010e+02
[136] -4.653843e+02 -9.966260e+01 -2.342686e+02  6.131400e+02  5.437831e+02
[141]  1.163759e+03 -3.602279e+02 -2.330141e+02  1.816029e+02  4.620475e+02
[146] -1.563333e+02 -9.049988e+01  2.980592e+02  1.583007e+02 -5.405435e+02
[151] -4.345100e+02 -2.262850e+02  7.099320e+02  5.437341e+02  2.500327e+02
[156] -2.808475e+01 -1.053728e+02  3.235875e+02 -4.429654e+02 -7.267447e+02
[161]  1.967229e+02  3.585873e+02 -3.737647e+02  4.723471e+02  2.016771e+02
[166]  5.890985e+02 -1.139407e+03  4.775878e+02  3.249846e+02 -2.190073e+02
[171] -7.443375e+02  4.077517e+02 -8.384774e+01  2.178708e+02  2.834072e+02
[176]  3.271277e+02  5.115315e+02 -3.510537e+02 -9.276523e+01  4.225447e+02
[181]  6.208844e+02  1.560188e+02 -6.478028e+01 -3.575554e+02  1.349416e+03
[186] -8.187161e+02 -2.873487e+02 -2.287618e+02 -7.404149e+02  7.848886e+02
[191] -3.954568e+02 -7.684468e+02 -1.893577e+02 -1.131051e+02  4.628189e+02
[196]  5.283408e+02  7.685769e+02  8.759307e+01 -6.300650e+02 -1.455912e+02
[201] -3.903437e+02 -8.834813e+02 -9.684689e+02 -8.582384e+02 -3.984895e+02
[206] -1.749766e+02  5.248759e+02  3.406808e+02  2.669995e+01 -4.779918e+02
[211]  3.094915e+02 -1.394421e+02 -8.864882e+01  5.356341e+02 -9.735456e+02
[216]  7.057131e+02  5.275957e+02 -1.511063e+02 -4.654692e+02  3.212737e+02
[221]  3.512970e+02  2.188433e+02  1.711722e+02  3.692600e+02 -3.821531e+02
[226] -5.090515e+02  1.828280e+02  2.474558e+02  1.707122e+02 -7.904289e+02
[231]  2.185227e+00  9.616593e+02 -8.304882e+02  5.448249e+02  1.327135e+03
[236]  5.914175e+01 -1.335112e+03  4.522697e+02 -3.379577e+03  1.428757e+03
[241]  2.391486e+02  6.021143e+02  3.725481e+02  5.272395e+02  4.007587e+02
[246] -2.393819e+02  5.496152e+02 -1.531555e+03 -1.348000e+03 -1.892949e+03
[251] -1.308407e+03  3.711359e+02  1.853660e+03  6.573426e+02  3.611924e+02
[256]  6.084509e+02  5.020492e+02 -8.523599e+02  5.949093e+02 -3.336415e+01
[261] -2.357101e+02  2.332169e+02 -6.910770e+02  6.463929e+02  7.331700e+02
[266] -1.989144e+03  1.875533e+03  1.747447e+02 -5.537421e+01 -4.533690e+02
[271] -3.842423e+02  9.376603e+02  4.816832e+02 -2.437853e+03 -6.530359e+02
[276]  8.664059e+01  5.770496e+02  4.495107e+02  1.265815e+02  9.708612e+02
[281]  9.416432e+02  6.454861e+02  4.056338e+01 -1.987265e+02 -2.107307e+03
[286] -7.123878e+02 -3.222680e+02  1.214500e+03  4.210129e+02  4.296327e+01
[291]  5.515583e+02 -2.017910e+03  7.379065e+02 -1.071835e+02 -1.146282e+02
[296] -9.580179e+01 -2.145587e+02  5.570525e+02 -4.152634e+02 -1.216878e+03
[301]  3.257690e+02 -3.304940e+03  6.946999e+02 -2.221130e+01  2.519095e+02
[306]  5.482436e+02  1.694075e+02  5.054563e+02 -3.791865e-02 -1.733648e+02
[311]  1.500267e+02  2.900246e+02 -2.476657e+01 -9.424374e+02 -8.761923e+01
[316]  8.130058e+02 -9.350384e+01  6.716697e+02  3.093363e+01 -2.136268e+03
[321]  1.459041e+02 -4.955556e+01  3.969771e+01 -2.210717e+02 -8.344575e+02
[326] -1.464870e+03 -6.727149e+01 -1.440194e+03  3.052321e+02  1.039981e+02
[331]  1.578326e+01  8.528166e+02 -6.023781e+02  8.973718e+02  4.837868e+02
[336]  5.788668e+02 -4.256644e+01 -2.001144e+01  4.393217e+02 -8.607459e+02
[341] -2.271823e+03  1.311481e+03  4.857369e+02 -2.870778e+02 -5.058675e+02
[346]  3.859803e+02  5.685632e+02  3.168231e+02  1.264277e+02  1.589492e+01
[351] -6.084816e+02 -7.088029e+02  3.955146e+02  3.604358e+02 -5.987423e+02
[356]  1.093632e+02 -9.194806e+02 -1.703859e+03 -1.388736e+03 -6.692028e+02
[361] -9.424814e+02  2.518577e+02 -4.187722e+00  6.165944e+02 -4.907684e+00
[366]  5.120736e-01 -3.191477e+02 -7.995134e+01  2.168426e+01  7.364285e+02
[371]  1.289494e+03  1.394570e+03  1.775672e+02 -5.460123e+02  1.182468e+03
[376] -1.058981e+03  1.347639e+03 -5.017082e+02 -6.725606e+02 -4.191570e+02
[381] -3.644357e+02  3.722062e+02  3.339692e+02  1.580054e+02 -9.813886e+01
[386] -1.835854e+03 -3.400440e+02  3.356695e+01  1.561257e+03  7.423095e+02
[391]  4.163702e+02  5.326623e+01  8.794563e+02 -2.105284e+02  2.554037e+02
[396]  1.032088e+03  5.493232e+02 -2.735709e+02  3.355016e+02 -8.132294e+02
[401] -2.042070e+02  5.335604e+02  4.861382e+02 -1.325424e+03  4.160361e+02
[406]  2.167781e+02 -1.595532e+03 -1.434754e+03  7.765817e+02  5.424896e+02
[411]  2.963024e+02 -8.329585e+02  9.723227e+02  8.522360e+02 -1.155845e+03
[416] -1.038581e+02  3.169473e+02  1.289166e+03  7.903067e+02 -9.583306e+02
[421] -7.117636e+02  4.944895e+02  8.158229e+02  1.410637e+02 -2.405275e+03
[426]  2.241828e+03 -9.694447e+02  4.258112e+02 -5.096736e+02 -1.596876e+02
[431]  6.471259e+02  8.536537e+02  9.604880e+02 -1.489535e+02  1.740700e+02
[436]  9.017514e+02  8.571139e+02  8.688850e+02  1.395133e+03  8.181980e+02
[441] -8.944033e+02  3.222456e+03 -3.511804e+02  6.561932e+02  4.278653e+02
[446]  4.865181e+02  1.032919e+03  1.904080e+03 -3.487265e+03  2.288547e+02
[451]  3.235525e+02 -7.656279e+02  9.548176e+00  1.103896e+02  3.928833e+01
[456]  5.080663e+02  1.069439e+02 -3.381635e+01  1.038880e+03  1.298233e+02
[461]  2.834813e+02 -6.722014e+01  8.453836e+02 -1.086434e+03 -6.591090e+01
[466]  1.644103e+02 -1.087635e+03 -1.218007e+02  7.072333e+02  1.454495e+03
[471]  4.943667e+02 -9.434966e+01  5.325716e+02 -1.946065e+03  1.157101e+03
[476]  8.092430e+02 -2.242046e+02 -5.334644e+03 -4.271690e+02  1.301632e+03
[481]  5.322536e+02 -9.368934e+02  8.959961e+02 -1.073603e+03  1.395001e+03
[486] -3.449438e+02 -6.932190e+01  9.233823e+02  5.769052e+02  2.815538e+02
[491]  5.728702e+02  6.253664e+02  5.541265e+02 -1.776640e+02 -1.287697e+03
[496]  1.264777e+03  6.898655e+02  1.086873e+03 -5.999267e+02 -3.154240e+03
[501]  7.281070e+02  1.862809e+03  5.809276e+02  7.741006e+02  1.565539e+03
[506]  2.798533e+02 -2.318901e+03  5.090777e+02 -8.641719e+02  9.054801e+02
[511] -5.104988e+01 -1.704701e+02  3.259073e+02  1.461158e+01 -7.545899e+01
[516]  6.703266e+02  7.466664e+02 -2.882277e+03  2.624908e+03  4.241838e+02
[521]  2.527210e+02  3.060668e+02  4.741113e+02  9.741708e+02  4.811258e+02
[526]  2.331294e+02 -6.403094e+02  1.202204e+02 -1.768983e+03  1.417322e+03
[531]  1.170967e+02  8.280731e+02  4.999849e+02 -3.263437e+02 -1.674329e+03
[536]  6.815213e+02 -5.860567e+02 -8.149747e+02 -6.071624e+02  1.104361e+03
[541]  1.058171e+02 -1.082636e+01  1.061093e+03  4.378186e+02 -1.664889e+01
[546] -1.545233e+03 -5.253680e+02 -7.313299e+02  2.775684e+02  2.254235e+02
[551]  9.025772e+02 -4.112427e+02  6.444017e+01 -1.410704e+03 -1.037925e+03
[556]  1.044071e+03 -2.527570e+02  9.480205e+02  7.154244e+02  1.000381e+03
[561]  3.242508e+02 -5.030482e+02  6.377183e+02  9.666666e+01 -1.100187e+03
[566]  2.922369e+02 -6.202441e+02 -1.549437e+03  2.126668e+03  2.090836e+02
[571]  1.018819e+03  1.186253e+03 -3.840128e+02  2.483740e+02 -1.957053e+02
[576] -5.939543e+01  3.219501e+02  2.445100e+02  6.505943e+02  1.527714e+02
[581]  2.212785e+02  9.699245e+01 -1.373256e+03  7.338067e+02  2.661148e+02
[586]  3.045275e+02  4.259434e+01 -1.210361e+03  1.122876e+02  2.581312e+01
[591]  1.762662e+02 -1.556148e+02  5.676017e+02  6.483222e+02  8.256401e+00
[596]  8.992129e+02 -2.618706e+03  8.375663e+02  3.046700e+02  3.590851e+02
[601]  5.579632e+02  4.059545e+02 -7.587075e+02 -1.177746e+03  9.005714e+02
[606]  1.083327e+02  7.705311e+02  3.479259e+02  2.381040e+02 -9.282413e+02
[611] -5.156819e+02 -2.107300e+02  3.725279e+02  2.460396e+02 -8.607810e+02
[616]  1.032533e+03 -9.246273e+02  2.132824e+03  1.611187e+02  6.326457e+02
[621]  5.999530e+02  3.239728e+02  6.422425e+02  1.084961e+03 -3.978300e+02
[626] -4.365553e+02 -2.931055e+03  1.703620e+03  4.540836e+02  5.416228e+02
[631]  9.037562e+02  2.268501e+02  1.805549e+02  2.202137e+01  2.583458e+02
[636] -3.454130e+02  5.205207e+00  8.984389e+02 -1.067874e+03 -4.413777e+02
[641] -2.098885e+03  1.635753e+03 -4.955828e+01  6.932879e+02  2.759365e+02
[646] -4.003513e+03  3.204611e+02  7.300354e+01  1.028199e+03  1.110674e+02
[651]  2.252434e+02  2.861661e+02 -1.723755e+02 -8.245870e+02  1.268322e+03
[656]  4.513397e+02  2.776592e+02 -1.846535e+03  1.731174e+03 -1.409509e+02
[661]  3.619413e+02  5.744842e+02  4.743343e+02  1.810784e+02  6.626898e+01
[666]  6.843578e+02 -3.013166e+03 -5.627797e+03 -2.659957e+03  1.661824e+03
[671] -5.717511e+01 -1.236017e+02 -4.343366e+02  4.057588e+02  3.564846e+02
[676]  1.226339e+02 -7.273028e+02  1.156895e+01  6.339418e+02  5.670390e+02
[681]  8.840298e+02  6.791554e+02 -1.052392e+03  7.722113e+02 -1.532320e+02
[686] -3.297784e+01 -1.116473e+02 -8.837177e+02  5.795412e+02  1.243874e+02
[691] -3.670301e+02 -2.912259e+03 -2.665485e+01 -2.327060e+03 -1.546317e+03
[696]  1.172590e+03 -7.421440e+02  1.264436e+03  5.448175e+02  1.025985e+03
[701]  1.662837e+02 -1.687331e+02  1.566771e+03  1.050726e+03 -7.780628e+01
[706]  9.435073e+01 -4.606934e+00  7.637658e+02 -1.932797e+03  7.578898e+02
[711]  4.388909e+02  3.794356e+01  2.428676e+02  2.011346e+02 -1.324570e+02
[716] -1.303193e+03  8.829519e+00  4.633785e+02 -9.777829e+01 -1.179380e+03
[721] -9.361207e+02 -2.493037e+03 -1.399341e+03 -2.379331e+03 -1.848254e+03
[726] -2.146511e+03 -5.623736e+01  4.067527e+02 -1.753871e+03 -1.465653e+02
[731]  6.067700e+02

$avar
[1] 829969.2

$xbar
[1] 4504.349

$aic
[1] 13.67292

$aicc
[1] 14.67683

$bic
[1] 13.77348

$se.phi
 [1] 0.001364571 0.001654490 0.001658689 0.001662209 0.001659617 0.001667593
 [7] 0.001690094 0.001693429 0.001691984 0.001664483 0.001662146 0.001659589
[13] 0.001661199 0.001662556 0.001396729

$se.theta
[1] 0
diff= artrans.wge(daily_bike_data$Total_Users, phi.tr = c(0,1))

aic5.wge(diff,type = 'bic') # bic picked p = 2 and q = 2
---------WORKING... PLEASE WAIT... 


Five Smallest Values of  bic 
    p    q        bic
    4    2   13.73395
    5    2   13.75966
    5    1   13.86287
    4    0   13.86295
    5    0   13.86810
seasonal_est = est.arma.wge(diff, p = 2, q = 2)
  
  
Coefficients of AR polynomial:  
0.3610 -0.0389 

                           AR Factor Table 
Factor                 Roots                Abs Recip    System Freq 
1-0.3610B+0.0389B^2    4.6361+-2.0477i      0.1973       0.0662
  
  
  
  
Coefficients of MA polynomial:  
-0.1184 0.8816 

                              MA FACTOR TABLE 
Factor                 Roots                Abs Recip    System Freq 
1+1.0000B             -1.0000               1.0000       0.5000
1-0.8816B              1.1343               0.8816       0.0000
  
  
# Checking the residuals to see if the model is appropriate 
plotts.wge(seasonal_est$res) # looks random
acf(seasonal_est$res) # all the acfs are within the bounds 
ljung.wge(seasonal_est$res) # greater than 0.05
Obs -0.5101017 0.5172595 -0.5272857 0.4848816 -0.512839 0.5486961 -0.4985657 0.513788 -0.5184669 0.5075917 -0.5210629 0.4958408 -0.5231253 0.5367516 -0.4706014 0.487064 -0.4970052 0.4831727 -0.5144025 0.505429 -0.4644367 0.5107547 -0.493999 0.4921956 
$test
[1] "Ljung-Box test"

$K
[1] 24

$chi.square
[1] 4569.486

$df
[1] 24

$pval
[1] 0
ljung.wge(seasonal_est$res, K =48) # greater than 0.05
Obs -0.5101017 0.5172595 -0.5272857 0.4848816 -0.512839 0.5486961 -0.4985657 0.513788 -0.5184669 0.5075917 -0.5210629 0.4958408 -0.5231253 0.5367516 -0.4706014 0.487064 -0.4970052 0.4831727 -0.5144025 0.505429 -0.4644367 0.5107547 -0.493999 0.4921956 -0.5112387 0.5169366 -0.4678605 0.518704 -0.4293692 0.4656065 -0.504262 0.4571246 -0.4830079 0.4933789 -0.4696618 0.4784432 -0.4553228 0.4490064 -0.468568 0.4628964 -0.4574949 0.4901441 -0.4715925 0.4778051 -0.5055763 0.4807059 -0.4659333 0.4842956 
$test
[1] "Ljung-Box test"

$K
[1] 48

$chi.square
[1] 8792.66

$df
[1] 48

$pval
[1] 0
seasonal_gen1 = gen.arima.wge(200, phi = seasonal_est$phi, theta = seasonal_est$theta, s = 2)
seasonal_gen2 = gen.arima.wge(200, phi = seasonal_est$phi, theta = seasonal_est$theta, s = 2)

seasonal_gen3 = gen.arima.wge(200, phi = seasonal_est$phi, theta = seasonal_est$theta, s = 2)
seasonal_gen4 = gen.arima.wge(200, phi = seasonal_est$phi, theta = seasonal_est$theta, s = 2)
plotts.sample.wge(seasonal_gen1)$xbar

[1] -2.034474
plotts.sample.wge(seasonal_gen2)$xbar

[1] -12.59861
plotts.sample.wge(seasonal_gen3)$xbar

[1] -10.57343
plotts.sample.wge(seasonal_gen4)$xbar # all these 4 random realizations have the same behavior as the original data so the model seems to be appropriate

[1] 1.590323
# Compare Spectral Densities: 
sims = 20
SpecDen = parzen.wge(daily_bike_data$Total_Users, plot = "FALSE")
plot(SpecDen$freq,SpecDen$pzgram, type = "l", lwd = 6)

for( i in 1: sims)
{
   SpecDen2 = parzen.wge(gen.aruma.wge(length(daily_bike_data$Total_Users), s = 2, phi = seasonal_est$phi, theta = seasonal_est$theta, plot ="FALSE"), plot = "FALSE")
   lines(SpecDen2$freq,SpecDen2$pzgram, lwd = 2, col = "red")
}

# Compare ACFs:
sims = 20
ACF = acf(daily_bike_data$Total_Users, plot = "FALSE")
plot(ACF$lag ,ACF$acf , type = "l", lwd = 6)

for( i in 1: sims)
{
   ACF2 = acf(gen.aruma.wge(length(daily_bike_data$Total_Users), s = 2, phi = seasonal_est$phi, theta = seasonal_est$theta, plot ="FALSE"), plot = "FALSE")
   lines(ACF2$lag ,ACF2$acf, lwd = 2, col = "red")
}

\[\phi(B)(1-B^2) = \theta(B)\] \[(1-0.1669356B-0.2335503B^2)(1-B^2) = (1- 0.005624988B- 0.823629022B^2)\]

# plot of Last N forecasts for short and long term horizon 
seasonal_stfore1 = fore.arima.wge(daily_bike_data$Total_Users, phi = seasonal_est$phi, theta = seasonal_est$theta, s = 2, n.ahead = 7, lastn = TRUE)
y.arma 364 761 251 44 -90 -647 -688 362 441 -159 143 259 -158 -217 -248 -521 650 1244 -107 -946 -557 435 999 -910 -1554 661 667 -71 403 264 25 190 182 -545 -85 707 -93 -107 8 141 -66 -157 441 226 202 660 812 -840 -1115 -528 -362 810 357 -456 162 941 -523 -551 688 -166 -190 392 -1339 -205 1528 19 -1510 86 1509 440 -86 -361 146 688 1047 373 -768 -1040 232 44 -838 89 631 -517 -468 732 -492 -740 691 567 1022 863 -1454 -307 1346 -1337 -686 1424 893 -861 -1186 1233 964 -2472 618 2634 -540 515 985 -2261 -153 2508 37 209 -201 -342 723 1254 -1244 -911 1100 -1768 -18 1975 281 -275 -352 470 -180 61 -77 -1455 448 549 -430 -103 452 1062 1230 -257 -1531 -168 704 185 -299 81 109 -660 -806 -124 986 1338 374 -406 -794 -73 -147 -918 185 1051 -126 54 431 160 -1124 -336 1352 -100 -1109 91 497 -45 484 412 314 -494 -657 517 867 137 -396 -713 924 16 -1414 -73 -589 744 841 -1250 -623 256 826 1196 839 -236 -1465 -761 -126 -757 -945 -499 219 555 984 816 -200 -810 85 456 -209 543 -692 -269 1292 -282 -1081 32 817 454 190 125 -642 -1085 188 905 356 -920 -541 1386 -280 -433 2022 372 -2353 -469 -2427 -327 3519 870 424 -89 -331 -631 213 -1133 -2230 -1355 -868 1548 3503 1502 -632 -283 72 -1104 -25 852 -486 28 -633 -187 1154 -1957 628 2615 -793 -890 -723 719 1295 -2410 -2284 1141 1538 1256 309 159 644 526 -292 -948 -2701 -1650 1228 2304 1397 -647 -293 -2146 -553 1880 113 77 -121 306 -293 -2028 -147 -2032 -416 3042 737 517 -94 -140 -48 -397 109 556 74 -1272 -741 1134 349 419 478 -2669 -1142 1575 610 128 -898 -1913 -199 -112 226 1573 279 799 -157 -254 813 327 -113 -455 197 -891 -3106 728 2915 -132 -877 120 780 430 186 -163 -970 -1146 664 1319 -743 -682 -451 -2057 -1455 306 408 985 1261 697 62 -705 -534 -58 417 1036 1730 1249 -673 -2145 173 -199 499 1037 -1604 -903 -195 624 1078 357 -213 -1991 -1186 1131 2362 1838 -264 -814 -52 -213 -399 1266 955 -748 -428 -929 -1204 952 1428 -982 -545 1029 -1661 -2302 1253 2393 747 -917 -15 1313 -1465 -1189 1088 1644 1285 -1286 -2330 -98 1590 974 -2488 627 1360 -924 229 -733 533 1583 1426 -347 -1264 342 1180 936 1014 345 -1934 1644 1514 -1683 201 77 778 2132 -3499 -3366 2186 106 140 1031 -239 102 582 -299 731 500 -315 24 400 -1291 -1272 749 -723 -509 1536 2051 734 -1090 -441 -2003 -126 2923 59 -6263 -3410 4606 2982 -607 37 -806 71 1352 -564 597 681 127 462 63 -610 -631 -1556 844 2313 857 -912 -4586 -1003 4581 2269 215 910 -510 -3935 -1056 901 697 1474 -234 -143 -493 -848 812 1595 -2728 782 3514 -1122 -640 57 493 681 4 -1138 -834 -1626 757 2391 244 339 -687 -2603 -153 1112 -920 -388 1553 1068 -679 551 556 -563 -1872 -1192 68 540 1129 1176 -419 -1196 -1401 -1535 1729 1618 695 1156 235 -477 -1468 -139 755 -1117 -195 157 -2132 1540 2507 182 1207 -731 -1269 -176 -307 420 619 475 45 -405 -437 -1711 189 1809 521 13 -1748 -987 758 584 240 464 821 -199 260 -2599 -1335 2457 845 759 207 -1712 -2327 864 1785 780 673 -347 -1573 -1540 -106 1054 1078 -661 392 -227 723 1549 -460 345 37 139 910 -676 -1845 -3260 722 3647 576 675 -260 -959 -369 297 -145 -318 1162 -526 -1777 -2250 794 2689 584 637 -4646 -2487 2882 2213 1178 -409 -461 -643 -1234 895 1586 -25 -2037 581 1400 -1032 642 635 -107 -249 493 -2985 -7830 -3363 5544 4890 281 -848 -740 121 579 -224 -371 957 1221 860 -267 -2758 -774 1351 203 184 -1029 -130 965 -353 -3209 -1236 -148 -1486 2810 1535 173 1364 408 -132 -1019 1043 1957 -505 -1231 -721 207 -1780 -412 2273 149 31 292 -485 -1825 -462 1771 682 -1429 -1644 -2379 -1836 -829 -774 -479 1101 2654 -773 -1299 1388 

seasonal_ltfore1 = fore.arima.wge(daily_bike_data$Total_Users, phi = seasonal_est$phi, theta = seasonal_est$theta, s = 2, n.ahead = 60, lastn = TRUE)
y.arma 364 761 251 44 -90 -647 -688 362 441 -159 143 259 -158 -217 -248 -521 650 1244 -107 -946 -557 435 999 -910 -1554 661 667 -71 403 264 25 190 182 -545 -85 707 -93 -107 8 141 -66 -157 441 226 202 660 812 -840 -1115 -528 -362 810 357 -456 162 941 -523 -551 688 -166 -190 392 -1339 -205 1528 19 -1510 86 1509 440 -86 -361 146 688 1047 373 -768 -1040 232 44 -838 89 631 -517 -468 732 -492 -740 691 567 1022 863 -1454 -307 1346 -1337 -686 1424 893 -861 -1186 1233 964 -2472 618 2634 -540 515 985 -2261 -153 2508 37 209 -201 -342 723 1254 -1244 -911 1100 -1768 -18 1975 281 -275 -352 470 -180 61 -77 -1455 448 549 -430 -103 452 1062 1230 -257 -1531 -168 704 185 -299 81 109 -660 -806 -124 986 1338 374 -406 -794 -73 -147 -918 185 1051 -126 54 431 160 -1124 -336 1352 -100 -1109 91 497 -45 484 412 314 -494 -657 517 867 137 -396 -713 924 16 -1414 -73 -589 744 841 -1250 -623 256 826 1196 839 -236 -1465 -761 -126 -757 -945 -499 219 555 984 816 -200 -810 85 456 -209 543 -692 -269 1292 -282 -1081 32 817 454 190 125 -642 -1085 188 905 356 -920 -541 1386 -280 -433 2022 372 -2353 -469 -2427 -327 3519 870 424 -89 -331 -631 213 -1133 -2230 -1355 -868 1548 3503 1502 -632 -283 72 -1104 -25 852 -486 28 -633 -187 1154 -1957 628 2615 -793 -890 -723 719 1295 -2410 -2284 1141 1538 1256 309 159 644 526 -292 -948 -2701 -1650 1228 2304 1397 -647 -293 -2146 -553 1880 113 77 -121 306 -293 -2028 -147 -2032 -416 3042 737 517 -94 -140 -48 -397 109 556 74 -1272 -741 1134 349 419 478 -2669 -1142 1575 610 128 -898 -1913 -199 -112 226 1573 279 799 -157 -254 813 327 -113 -455 197 -891 -3106 728 2915 -132 -877 120 780 430 186 -163 -970 -1146 664 1319 -743 -682 -451 -2057 -1455 306 408 985 1261 697 62 -705 -534 -58 417 1036 1730 1249 -673 -2145 173 -199 499 1037 -1604 -903 -195 624 1078 357 -213 -1991 -1186 1131 2362 1838 -264 -814 -52 -213 -399 1266 955 -748 -428 -929 -1204 952 1428 -982 -545 1029 -1661 -2302 1253 2393 747 -917 -15 1313 -1465 -1189 1088 1644 1285 -1286 -2330 -98 1590 974 -2488 627 1360 -924 229 -733 533 1583 1426 -347 -1264 342 1180 936 1014 345 -1934 1644 1514 -1683 201 77 778 2132 -3499 -3366 2186 106 140 1031 -239 102 582 -299 731 500 -315 24 400 -1291 -1272 749 -723 -509 1536 2051 734 -1090 -441 -2003 -126 2923 59 -6263 -3410 4606 2982 -607 37 -806 71 1352 -564 597 681 127 462 63 -610 -631 -1556 844 2313 857 -912 -4586 -1003 4581 2269 215 910 -510 -3935 -1056 901 697 1474 -234 -143 -493 -848 812 1595 -2728 782 3514 -1122 -640 57 493 681 4 -1138 -834 -1626 757 2391 244 339 -687 -2603 -153 1112 -920 -388 1553 1068 -679 551 556 -563 -1872 -1192 68 540 1129 1176 -419 -1196 -1401 -1535 1729 1618 695 1156 235 -477 -1468 -139 755 -1117 -195 157 -2132 1540 2507 182 1207 -731 -1269 -176 -307 420 619 475 45 -405 -437 -1711 189 1809 521 13 -1748 -987 758 584 240 464 821 -199 260 -2599 -1335 2457 845 759 207 -1712 -2327 864 1785 780 673 -347 -1573 -1540 -106 1054 1078 -661 392 -227 723 1549 -460 345 37 139 910 -676 -1845 -3260 722 3647 576 675 -260 -959 -369 297 -145 -318 1162 -526 -1777 -2250 794 2689 584 637 -4646 -2487 2882 2213 1178 -409 -461 -643 -1234 895 1586 -25 -2037 581 1400 -1032 642 635 -107 -249 493 -2985 -7830 -3363 5544 4890 281 -848 -740 121 579 -224 -371 957 1221 860 -267 -2758 -774 1351 203 184 -1029 -130 965 -353 -3209 -1236 -148 -1486 2810 1535 173 1364 408 -132 -1019 1043 1957 -505 -1231 -721 207 -1780 -412 2273 149 31 292 -485 -1825 -462 1771 682 -1429 -1644 -2379 -1836 -829 -774 -479 1101 2654 -773 -1299 1388 

t = 1:length(daily_bike_data$Total_Users)
plot(t[720:731],daily_bike_data$Total_Users[720:731], type = 'l', xlab = "Time", ylab = "Total Users")
points(t[725:731], seasonal_stfore1$f, type="l", lwd=2, lty = 2, col = 'blue')

plot(t[670:731],daily_bike_data$Total_Users[670:731], type = 'l', xlab = "Time", ylab = "Total Users")
points(t[672:731], seasonal_ltfore1$f, type="l", lwd=2, lty = 2,col = 'blue')

# ASE 
seasonal_stASE = mean((daily_bike_data$Total_Users[725:731]-seasonal_stfore1$f)^2)
seasonal_ltASE = mean((daily_bike_data$Total_Users[672:731]-seasonal_ltfore1$f)^2)
seasonal_stASE
[1] 4368663
seasonal_ltASE
[1] 4042180
# Rolling Window RMSE
# RW-RMSE commented out due to obscene amount of unsupressable output
# seasonal_strwRMSE = roll.win.rmse.wge(daily_bike_data$Total_Users, horizon = 7, phi = seasonal_est$phi, theta = seasonal_est$theta, s = 2)$rwRMSE
# seasonal_ltrwRMSE = roll.win.rmse.wge(daily_bike_data$Total_Users, horizon = 60, phi = seasonal_est$phi, theta = seasonal_est$theta, s = 2)$rwRMSE

seasonal_strwRMSE = 1253.353
seasonal_ltrwRMSE = 1504.218

seasonal_stfore2 = fore.arima.wge(daily_bike_data$Total_Users, phi = seasonal_est$phi, theta = seasonal_est$theta, s = 2, n.ahead = 7)
y.arma 364 761 251 44 -90 -647 -688 362 441 -159 143 259 -158 -217 -248 -521 650 1244 -107 -946 -557 435 999 -910 -1554 661 667 -71 403 264 25 190 182 -545 -85 707 -93 -107 8 141 -66 -157 441 226 202 660 812 -840 -1115 -528 -362 810 357 -456 162 941 -523 -551 688 -166 -190 392 -1339 -205 1528 19 -1510 86 1509 440 -86 -361 146 688 1047 373 -768 -1040 232 44 -838 89 631 -517 -468 732 -492 -740 691 567 1022 863 -1454 -307 1346 -1337 -686 1424 893 -861 -1186 1233 964 -2472 618 2634 -540 515 985 -2261 -153 2508 37 209 -201 -342 723 1254 -1244 -911 1100 -1768 -18 1975 281 -275 -352 470 -180 61 -77 -1455 448 549 -430 -103 452 1062 1230 -257 -1531 -168 704 185 -299 81 109 -660 -806 -124 986 1338 374 -406 -794 -73 -147 -918 185 1051 -126 54 431 160 -1124 -336 1352 -100 -1109 91 497 -45 484 412 314 -494 -657 517 867 137 -396 -713 924 16 -1414 -73 -589 744 841 -1250 -623 256 826 1196 839 -236 -1465 -761 -126 -757 -945 -499 219 555 984 816 -200 -810 85 456 -209 543 -692 -269 1292 -282 -1081 32 817 454 190 125 -642 -1085 188 905 356 -920 -541 1386 -280 -433 2022 372 -2353 -469 -2427 -327 3519 870 424 -89 -331 -631 213 -1133 -2230 -1355 -868 1548 3503 1502 -632 -283 72 -1104 -25 852 -486 28 -633 -187 1154 -1957 628 2615 -793 -890 -723 719 1295 -2410 -2284 1141 1538 1256 309 159 644 526 -292 -948 -2701 -1650 1228 2304 1397 -647 -293 -2146 -553 1880 113 77 -121 306 -293 -2028 -147 -2032 -416 3042 737 517 -94 -140 -48 -397 109 556 74 -1272 -741 1134 349 419 478 -2669 -1142 1575 610 128 -898 -1913 -199 -112 226 1573 279 799 -157 -254 813 327 -113 -455 197 -891 -3106 728 2915 -132 -877 120 780 430 186 -163 -970 -1146 664 1319 -743 -682 -451 -2057 -1455 306 408 985 1261 697 62 -705 -534 -58 417 1036 1730 1249 -673 -2145 173 -199 499 1037 -1604 -903 -195 624 1078 357 -213 -1991 -1186 1131 2362 1838 -264 -814 -52 -213 -399 1266 955 -748 -428 -929 -1204 952 1428 -982 -545 1029 -1661 -2302 1253 2393 747 -917 -15 1313 -1465 -1189 1088 1644 1285 -1286 -2330 -98 1590 974 -2488 627 1360 -924 229 -733 533 1583 1426 -347 -1264 342 1180 936 1014 345 -1934 1644 1514 -1683 201 77 778 2132 -3499 -3366 2186 106 140 1031 -239 102 582 -299 731 500 -315 24 400 -1291 -1272 749 -723 -509 1536 2051 734 -1090 -441 -2003 -126 2923 59 -6263 -3410 4606 2982 -607 37 -806 71 1352 -564 597 681 127 462 63 -610 -631 -1556 844 2313 857 -912 -4586 -1003 4581 2269 215 910 -510 -3935 -1056 901 697 1474 -234 -143 -493 -848 812 1595 -2728 782 3514 -1122 -640 57 493 681 4 -1138 -834 -1626 757 2391 244 339 -687 -2603 -153 1112 -920 -388 1553 1068 -679 551 556 -563 -1872 -1192 68 540 1129 1176 -419 -1196 -1401 -1535 1729 1618 695 1156 235 -477 -1468 -139 755 -1117 -195 157 -2132 1540 2507 182 1207 -731 -1269 -176 -307 420 619 475 45 -405 -437 -1711 189 1809 521 13 -1748 -987 758 584 240 464 821 -199 260 -2599 -1335 2457 845 759 207 -1712 -2327 864 1785 780 673 -347 -1573 -1540 -106 1054 1078 -661 392 -227 723 1549 -460 345 37 139 910 -676 -1845 -3260 722 3647 576 675 -260 -959 -369 297 -145 -318 1162 -526 -1777 -2250 794 2689 584 637 -4646 -2487 2882 2213 1178 -409 -461 -643 -1234 895 1586 -25 -2037 581 1400 -1032 642 635 -107 -249 493 -2985 -7830 -3363 5544 4890 281 -848 -740 121 579 -224 -371 957 1221 860 -267 -2758 -774 1351 203 184 -1029 -130 965 -353 -3209 -1236 -148 -1486 2810 1535 173 1364 408 -132 -1019 1043 1957 -505 -1231 -721 207 -1780 -412 2273 149 31 292 -485 -1825 -462 1771 682 -1429 -1644 -2379 -1836 -829 -774 -479 1101 2654 -773 -1299 1388 

seasonal_ltfore2 = fore.arima.wge(daily_bike_data$Total_Users, phi = seasonal_est$phi, theta = seasonal_est$theta, s = 2, n.ahead = 60)
y.arma 364 761 251 44 -90 -647 -688 362 441 -159 143 259 -158 -217 -248 -521 650 1244 -107 -946 -557 435 999 -910 -1554 661 667 -71 403 264 25 190 182 -545 -85 707 -93 -107 8 141 -66 -157 441 226 202 660 812 -840 -1115 -528 -362 810 357 -456 162 941 -523 -551 688 -166 -190 392 -1339 -205 1528 19 -1510 86 1509 440 -86 -361 146 688 1047 373 -768 -1040 232 44 -838 89 631 -517 -468 732 -492 -740 691 567 1022 863 -1454 -307 1346 -1337 -686 1424 893 -861 -1186 1233 964 -2472 618 2634 -540 515 985 -2261 -153 2508 37 209 -201 -342 723 1254 -1244 -911 1100 -1768 -18 1975 281 -275 -352 470 -180 61 -77 -1455 448 549 -430 -103 452 1062 1230 -257 -1531 -168 704 185 -299 81 109 -660 -806 -124 986 1338 374 -406 -794 -73 -147 -918 185 1051 -126 54 431 160 -1124 -336 1352 -100 -1109 91 497 -45 484 412 314 -494 -657 517 867 137 -396 -713 924 16 -1414 -73 -589 744 841 -1250 -623 256 826 1196 839 -236 -1465 -761 -126 -757 -945 -499 219 555 984 816 -200 -810 85 456 -209 543 -692 -269 1292 -282 -1081 32 817 454 190 125 -642 -1085 188 905 356 -920 -541 1386 -280 -433 2022 372 -2353 -469 -2427 -327 3519 870 424 -89 -331 -631 213 -1133 -2230 -1355 -868 1548 3503 1502 -632 -283 72 -1104 -25 852 -486 28 -633 -187 1154 -1957 628 2615 -793 -890 -723 719 1295 -2410 -2284 1141 1538 1256 309 159 644 526 -292 -948 -2701 -1650 1228 2304 1397 -647 -293 -2146 -553 1880 113 77 -121 306 -293 -2028 -147 -2032 -416 3042 737 517 -94 -140 -48 -397 109 556 74 -1272 -741 1134 349 419 478 -2669 -1142 1575 610 128 -898 -1913 -199 -112 226 1573 279 799 -157 -254 813 327 -113 -455 197 -891 -3106 728 2915 -132 -877 120 780 430 186 -163 -970 -1146 664 1319 -743 -682 -451 -2057 -1455 306 408 985 1261 697 62 -705 -534 -58 417 1036 1730 1249 -673 -2145 173 -199 499 1037 -1604 -903 -195 624 1078 357 -213 -1991 -1186 1131 2362 1838 -264 -814 -52 -213 -399 1266 955 -748 -428 -929 -1204 952 1428 -982 -545 1029 -1661 -2302 1253 2393 747 -917 -15 1313 -1465 -1189 1088 1644 1285 -1286 -2330 -98 1590 974 -2488 627 1360 -924 229 -733 533 1583 1426 -347 -1264 342 1180 936 1014 345 -1934 1644 1514 -1683 201 77 778 2132 -3499 -3366 2186 106 140 1031 -239 102 582 -299 731 500 -315 24 400 -1291 -1272 749 -723 -509 1536 2051 734 -1090 -441 -2003 -126 2923 59 -6263 -3410 4606 2982 -607 37 -806 71 1352 -564 597 681 127 462 63 -610 -631 -1556 844 2313 857 -912 -4586 -1003 4581 2269 215 910 -510 -3935 -1056 901 697 1474 -234 -143 -493 -848 812 1595 -2728 782 3514 -1122 -640 57 493 681 4 -1138 -834 -1626 757 2391 244 339 -687 -2603 -153 1112 -920 -388 1553 1068 -679 551 556 -563 -1872 -1192 68 540 1129 1176 -419 -1196 -1401 -1535 1729 1618 695 1156 235 -477 -1468 -139 755 -1117 -195 157 -2132 1540 2507 182 1207 -731 -1269 -176 -307 420 619 475 45 -405 -437 -1711 189 1809 521 13 -1748 -987 758 584 240 464 821 -199 260 -2599 -1335 2457 845 759 207 -1712 -2327 864 1785 780 673 -347 -1573 -1540 -106 1054 1078 -661 392 -227 723 1549 -460 345 37 139 910 -676 -1845 -3260 722 3647 576 675 -260 -959 -369 297 -145 -318 1162 -526 -1777 -2250 794 2689 584 637 -4646 -2487 2882 2213 1178 -409 -461 -643 -1234 895 1586 -25 -2037 581 1400 -1032 642 635 -107 -249 493 -2985 -7830 -3363 5544 4890 281 -848 -740 121 579 -224 -371 957 1221 860 -267 -2758 -774 1351 203 184 -1029 -130 965 -353 -3209 -1236 -148 -1486 2810 1535 173 1364 408 -132 -1019 1043 1957 -505 -1231 -721 207 -1780 -412 2273 149 31 292 -485 -1825 -462 1771 682 -1429 -1644 -2379 -1836 -829 -774 -479 1101 2654 -773 -1299 1388 

# Plots of the short and Long Term Forecasts  
t = 1:800
plot(t[670:731],daily_bike_data$Total_Users[670:731], type = 'l', main = "Short Term Forecast", xlim = c(670,795), xlab = "Time", ylab = "Total Users")
points(t[732:738],seasonal_stfore2$f, type = 'l', col = 'blue')
points(t[732:738],seasonal_stfore2$ll, type = 'l',lwd=2, lty = 2, col = 'red')
points(t[732:738],seasonal_stfore2$ul, type = 'l',lwd=2, lty = 2, col = 'red')

plot(t[670:731],daily_bike_data$Total_Users[670:731], type = 'l', main = "Long Term Forecast", xlim = c(670,795), xlab = "Time", ylab = "Total Users")
points(t[732:791],seasonal_ltfore2$f, type = 'l', col = 'blue')
points(t[732:791],seasonal_ltfore2$ll, type = 'l',lwd=2, lty = 2, col = 'red')
points(t[732:791],seasonal_ltfore2$ul, type = 'l',lwd=2, lty = 2, col = 'red')

Non-Seasonal ARIMA Model

# Check if the data could be stationary
adf.test(daily_bike_data$Total_Users) # p-value 0.7, data not likely stationary

    Augmented Dickey-Fuller Test

data:  daily_bike_data$Total_Users
Dickey-Fuller = -1.6351, Lag order = 9, p-value = 0.7327
alternative hypothesis: stationary
# Difference of 1
total_d1 = artrans.wge(daily_bike_data$Total_Users, phi.tr = c(1))

# Model the residuals
aic5.wge(total_d1, type = 'bic') # bic and aic pick p = 1 and q = 1
---------WORKING... PLEASE WAIT... 


Five Smallest Values of  bic 
    p    q        bic
    1    1   13.68620
    0    2   13.69304
    2    1   13.69387
    1    2   13.69455
    3    1   13.69895
est = est.arma.wge(total_d1, p = 1, q = 1)
  
  
Coefficients of AR polynomial:  
0.3591 

                           AR Factor Table 
Factor                 Roots                Abs Recip    System Freq 
1-0.3591B              2.7848               0.3591       0.0000
  
  
  
  
Coefficients of MA polynomial:  
0.8903 

                              MA FACTOR TABLE 
Factor                 Roots                Abs Recip    System Freq 
1-0.8903B              1.1232               0.8903       0.0000
  
  
plotts.sample.wge(est$res, arlimits = TRUE)$xbar # appears to be white noise

[1] -1.083662
ljung.wge(est$res) # FTR
Obs 0.01106095 -0.01226944 -0.04955624 -0.04487202 -0.01066982 0.09196727 0.01823558 0.01832273 -0.0280327 0.008467002 -0.03320097 -0.01648752 -0.04075165 0.07603663 0.06427303 -0.0260844 0.005665222 -0.03280932 -0.03354154 0.01664482 0.07094442 0.0320215 -0.0004189083 -0.006244533 
$test
[1] "Ljung-Box test"

$K
[1] 24

$chi.square
[1] 27.5971

$df
[1] 24

$pval
[1] 0.2773996
ljung.wge(est$res, K = 48) # Reject. Mixed evidence.
Obs 0.01106095 -0.01226944 -0.04955624 -0.04487202 -0.01066982 0.09196727 0.01823558 0.01832273 -0.0280327 0.008467002 -0.03320097 -0.01648752 -0.04075165 0.07603663 0.06427303 -0.0260844 0.005665222 -0.03280932 -0.03354154 0.01664482 0.07094442 0.0320215 -0.0004189083 -0.006244533 -0.03861073 0.04916367 0.04703067 0.06298134 0.1247115 -0.04713436 -0.03531116 -0.06145962 0.01091111 0.0188672 0.03310521 -0.006645202 0.05320788 -0.06603257 0.0242021 -0.03893266 0.05058829 0.02162102 0.02101459 -0.007107952 -0.05170156 0.002141853 0.02721173 0.01534891 
$test
[1] "Ljung-Box test"

$K
[1] 48

$chi.square
[1] 66.66045

$df
[1] 48

$pval
[1] 0.03852057
# Generate realizations for comparison
arima_gen1 = gen.arima.wge(200, phi = est$phi, theta = est$theta, d = 1)

arima_gen2 = gen.arima.wge(200, phi = est$phi, theta = est$theta, d = 1)

arima_gen3 = gen.arima.wge(200, phi = est$phi, theta = est$theta, d = 1)

arima_gen4 = gen.arima.wge(200, phi = est$phi, theta = est$theta, d = 1)

plotts.sample.wge(arima_gen1)$xbar

[1] 11.9336
plotts.sample.wge(arima_gen2)$xbar

[1] -3.912301
plotts.sample.wge(arima_gen3)$xbar

[1] -11.69953
plotts.sample.wge(arima_gen4)$xbar # Similar behavior to original data

[1] -0.9305175
# Compare Spectral Densities: 
sims = 20
SpecDen = parzen.wge(daily_bike_data$Total_Users, plot = "FALSE")
plot(SpecDen$freq,SpecDen$pzgram, type = "l", lwd = 6)

for( i in 1: sims)
{
   SpecDen2 = parzen.wge(gen.aruma.wge(length(daily_bike_data$Total_Users), d = 1, phi = est$phi, theta = est$theta, plot ="FALSE"), plot = "FALSE")
   lines(SpecDen2$freq,SpecDen2$pzgram, lwd = 2, col = "red")
}

# Compare ACFs:
sims = 20
ACF = acf(daily_bike_data$Total_Users, plot = "FALSE")
plot(ACF$lag ,ACF$acf , type = "l", lwd = 6)

for( i in 1: sims)
{
   ACF2 = acf(gen.aruma.wge(length(daily_bike_data$Total_Users), d = 1, phi = est$phi, theta = est$theta, plot ="FALSE"), plot = "FALSE")
   lines(ACF2$lag ,ACF2$acf, lwd = 2, col = "red")
}

# Short and Long Term Forecasts
fore_arima_st = fore.arima.wge(daily_bike_data$Total_Users, phi = est$phi, theta = est$theta, d = 1, n.ahead = 7, lastn = TRUE)
y.arma -184 548 213 38 6 -96 -551 -137 499 -58 -101 244 15 -173 -44 -204 -317 967 277 -384 -562 5 430 569 -1479 -75 736 -69 -2 405 -141 166 24 158 -703 618 89 -182 75 -67 208 -274 117 324 -98 300 360 452 -1292 177 -705 343 467 -110 -346 508 433 -956 405 283 -449 259 133 -1472 1267 261 -242 -1268 1354 155 285 -371 10 136 552 495 -122 -646 -394 626 -582 -256 345 286 -803 335 397 -889 149 542 25 997 -134 -1320 1013 333 -1670 984 440 453 -1314 128 1105 -141 -2331 2949 -315 -225 740 245 -2506 2353 155 -118 327 -528 186 537 717 -1961 1050 50 -1818 1800 175 106 -381 29 441 -621 682 -759 -696 1144 -595 165 -268 720 342 888 -1145 -386 218 486 -301 2 79 30 -690 -116 -8 994 344 30 -436 -358 285 -432 -486 671 380 -506 560 -129 289 -1413 1077 275 -375 -734 825 -328 283 201 211 103 -597 -60 577 290 -153 -243 -470 1394 -1378 -36 -37 -552 1296 -455 -795 172 84 742 454 385 -621 -844 83 -209 -548 -397 -102 321 234 750 66 -266 -544 629 -173 -36 579 -1271 1002 290 -572 -509 541 276 178 12 113 -755 -330 518 387 -31 -889 348 1038 -1318 885 1137 -765 -1588 1119 -3546 3219 300 570 -146 57 -388 -243 456 -1589 -641 -714 -154 1702 1801 -299 -333 50 22 -1126 1101 -249 -237 265 -898 711 443 -2400 3028 -413 -380 -510 -213 932 363 -2773 489 652 886 370 -61 220 424 102 -394 -554 -2147 497 731 1573 -176 -471 178 -2324 1771 109 4 73 -194 500 -793 -1235 1088 -3120 2704 338 399 118 -212 72 -120 -277 386 170 -96 -1176 435 699 -350 769 -291 -2378 1236 339 271 -143 -755 -1158 959 -1071 1297 276 3 796 -953 699 114 213 -326 -129 326 -1217 -1889 2617 298 -430 -447 567 213 217 -31 -132 -838 -308 972 347 -1090 408 -859 -1198 -257 563 -155 1140 121 576 -514 -191 -343 285 132 904 826 423 -1096 -1049 1222 -1421 1920 -883 -721 -182 -13 637 441 -84 -129 -1862 676 455 1907 -69 -195 -619 567 -780 381 885 70 -818 390 -1319 115 837 591 -1573 1028 1 -1662 -640 1893 500 247 -1164 1149 164 -1629 440 648 996 289 -1575 -755 657 933 41 -2529 3156 -1796 872 -643 -90 623 960 466 -813 -451 793 387 549 465 -120 -1814 3458 -1944 261 -60 137 641 1491 -4990 1624 562 -456 596 435 -674 776 -194 -105 836 -336 21 3 397 -1688 416 333 -1056 547 989 1062 -328 -762 321 -2324 2198 725 -666 -5597 2187 2419 563 -1170 1207 -2013 2084 -732 168 429 252 -125 587 -524 -86 -545 -1011 1855 458 399 -1311 -3275 2272 2309 -40 255 655 -1165 -2770 1714 -813 1510 -36 -198 55 -548 -300 1112 483 -3211 3993 -479 -643 3 54 439 242 -238 -900 66 -1692 2449 -58 302 37 -724 -1879 1726 -614 -306 -82 1635 -567 -112 663 -107 -456 -1416 224 -156 696 433 743 -1162 -34 -1367 -168 1897 -279 974 182 53 -530 -938 799 -44 -1073 878 -721 -1411 2951 -444 626 581 -1312 43 -219 -88 508 111 364 -319 -86 -351 -1360 1549 260 261 -248 -1500 513 245 339 -99 563 258 -457 717 -3316 1981 476 369 390 -183 -1529 -798 1662 123 657 16 -363 -1210 -330 224 830 248 -909 1301 -1528 2251 -702 242 103 -66 205 705 -1381 -464 -2796 3518 129 447 228 -488 -471 102 195 -340 22 1140 -1666 -111 -2139 2933 -244 828 -191 -4455 1968 914 1299 -121 -288 -173 -470 -764 1659 -73 48 -2085 2666 -1266 234 408 227 -334 85 408 -3393 -4437 1074 4470 420 -139 -709 -31 152 427 -651 280 677 544 316 -583 -2175 1401 -50 253 -69 -960 830 135 -488 -2721 1485 -1633 147 2663 -1128 1301 63 345 -477 -542 1585 372 -877 -354 -367 574 -2354 1942 331 -182 213 79 -564 -1261 799 972 -290 -1139 -505 -1874 38 -867 93 -572 1673 981 -1754 455 933 

fore_arima_lt = fore.arima.wge(daily_bike_data$Total_Users, phi = est$phi, theta = est$theta, d = 1, n.ahead = 60, lastn = TRUE)
y.arma -184 548 213 38 6 -96 -551 -137 499 -58 -101 244 15 -173 -44 -204 -317 967 277 -384 -562 5 430 569 -1479 -75 736 -69 -2 405 -141 166 24 158 -703 618 89 -182 75 -67 208 -274 117 324 -98 300 360 452 -1292 177 -705 343 467 -110 -346 508 433 -956 405 283 -449 259 133 -1472 1267 261 -242 -1268 1354 155 285 -371 10 136 552 495 -122 -646 -394 626 -582 -256 345 286 -803 335 397 -889 149 542 25 997 -134 -1320 1013 333 -1670 984 440 453 -1314 128 1105 -141 -2331 2949 -315 -225 740 245 -2506 2353 155 -118 327 -528 186 537 717 -1961 1050 50 -1818 1800 175 106 -381 29 441 -621 682 -759 -696 1144 -595 165 -268 720 342 888 -1145 -386 218 486 -301 2 79 30 -690 -116 -8 994 344 30 -436 -358 285 -432 -486 671 380 -506 560 -129 289 -1413 1077 275 -375 -734 825 -328 283 201 211 103 -597 -60 577 290 -153 -243 -470 1394 -1378 -36 -37 -552 1296 -455 -795 172 84 742 454 385 -621 -844 83 -209 -548 -397 -102 321 234 750 66 -266 -544 629 -173 -36 579 -1271 1002 290 -572 -509 541 276 178 12 113 -755 -330 518 387 -31 -889 348 1038 -1318 885 1137 -765 -1588 1119 -3546 3219 300 570 -146 57 -388 -243 456 -1589 -641 -714 -154 1702 1801 -299 -333 50 22 -1126 1101 -249 -237 265 -898 711 443 -2400 3028 -413 -380 -510 -213 932 363 -2773 489 652 886 370 -61 220 424 102 -394 -554 -2147 497 731 1573 -176 -471 178 -2324 1771 109 4 73 -194 500 -793 -1235 1088 -3120 2704 338 399 118 -212 72 -120 -277 386 170 -96 -1176 435 699 -350 769 -291 -2378 1236 339 271 -143 -755 -1158 959 -1071 1297 276 3 796 -953 699 114 213 -326 -129 326 -1217 -1889 2617 298 -430 -447 567 213 217 -31 -132 -838 -308 972 347 -1090 408 -859 -1198 -257 563 -155 1140 121 576 -514 -191 -343 285 132 904 826 423 -1096 -1049 1222 -1421 1920 -883 -721 -182 -13 637 441 -84 -129 -1862 676 455 1907 -69 -195 -619 567 -780 381 885 70 -818 390 -1319 115 837 591 -1573 1028 1 -1662 -640 1893 500 247 -1164 1149 164 -1629 440 648 996 289 -1575 -755 657 933 41 -2529 3156 -1796 872 -643 -90 623 960 466 -813 -451 793 387 549 465 -120 -1814 3458 -1944 261 -60 137 641 1491 -4990 1624 562 -456 596 435 -674 776 -194 -105 836 -336 21 3 397 -1688 416 333 -1056 547 989 1062 -328 -762 321 -2324 2198 725 -666 -5597 2187 2419 563 -1170 1207 -2013 2084 -732 168 429 252 -125 587 -524 -86 -545 -1011 1855 458 399 -1311 -3275 2272 2309 -40 255 655 -1165 -2770 1714 -813 1510 -36 -198 55 -548 -300 1112 483 -3211 3993 -479 -643 3 54 439 242 -238 -900 66 -1692 2449 -58 302 37 -724 -1879 1726 -614 -306 -82 1635 -567 -112 663 -107 -456 -1416 224 -156 696 433 743 -1162 -34 -1367 -168 1897 -279 974 182 53 -530 -938 799 -44 -1073 878 -721 -1411 2951 -444 626 581 -1312 43 -219 -88 508 111 364 -319 -86 -351 -1360 1549 260 261 -248 -1500 513 245 339 -99 563 258 -457 717 -3316 1981 476 369 390 -183 -1529 -798 1662 123 657 16 -363 -1210 -330 224 830 248 -909 1301 -1528 2251 -702 242 103 -66 205 705 -1381 -464 -2796 3518 129 447 228 -488 -471 102 195 -340 22 1140 -1666 -111 -2139 2933 -244 828 -191 -4455 1968 914 1299 -121 -288 -173 -470 -764 1659 -73 48 -2085 2666 -1266 234 408 227 -334 85 408 -3393 -4437 1074 4470 420 -139 -709 -31 152 427 -651 280 677 544 316 -583 -2175 1401 -50 253 -69 -960 830 135 -488 -2721 1485 -1633 147 2663 -1128 1301 63 345 -477 -542 1585 372 -877 -354 -367 574 -2354 1942 331 -182 213 79 -564 -1261 799 972 -290 -1139 -505 -1874 38 -867 93 -572 1673 981 -1754 455 933 

# Plot Forecasts and CI
t = 1:length(daily_bike_data$Total_Users)
plot(t[720:731],daily_bike_data$Total_Users[720:731], type = 'l', xlab = "Time", ylab = "Total Users")
points(t[725:731], fore_arima_st$f, type="l", lwd=2, lty = 2, col = 'blue')

plot(t[670:731],daily_bike_data$Total_Users[670:731], type = 'l', xlab = "Time", ylab = "Total Users")
points(t[672:731], fore_arima_lt$f, type="l", lwd=2, lty = 2,col = 'blue')

# ASE 
arima_stASE = mean((daily_bike_data$Total_Users[725:731]-fore_arima_st$f)^2)
arima_ltASE = mean((daily_bike_data$Total_Users[672:731]-fore_arima_lt$f)^2)
arima_stASE
[1] 3230238
arima_ltASE
[1] 3930824
# Rolling Window RMSE
# RW-RMSE commented out due to obscene amount of unsupressable output
# arima_strwRMSE = roll.win.rmse.wge(daily_bike_data$Total_Users, horizon = 7, phi = est$phi, theta = est$theta, d = 1)$rwRMSE
# arima_ltrwRMSE = roll.win.rmse.wge(daily_bike_data$Total_Users, horizon = 60, phi = est$phi, theta = est$theta, d = 1)$rwRMSE
arima_strwRMSE = 1237.393
arima_ltrwRMSE = 1503.197

# Plots of the short and Long Term Forecasts  
fore_arima_st_2 = fore.arima.wge(daily_bike_data$Total_Users, phi = est$phi, theta = est$theta, d = 1, n.ahead = 7)
y.arma -184 548 213 38 6 -96 -551 -137 499 -58 -101 244 15 -173 -44 -204 -317 967 277 -384 -562 5 430 569 -1479 -75 736 -69 -2 405 -141 166 24 158 -703 618 89 -182 75 -67 208 -274 117 324 -98 300 360 452 -1292 177 -705 343 467 -110 -346 508 433 -956 405 283 -449 259 133 -1472 1267 261 -242 -1268 1354 155 285 -371 10 136 552 495 -122 -646 -394 626 -582 -256 345 286 -803 335 397 -889 149 542 25 997 -134 -1320 1013 333 -1670 984 440 453 -1314 128 1105 -141 -2331 2949 -315 -225 740 245 -2506 2353 155 -118 327 -528 186 537 717 -1961 1050 50 -1818 1800 175 106 -381 29 441 -621 682 -759 -696 1144 -595 165 -268 720 342 888 -1145 -386 218 486 -301 2 79 30 -690 -116 -8 994 344 30 -436 -358 285 -432 -486 671 380 -506 560 -129 289 -1413 1077 275 -375 -734 825 -328 283 201 211 103 -597 -60 577 290 -153 -243 -470 1394 -1378 -36 -37 -552 1296 -455 -795 172 84 742 454 385 -621 -844 83 -209 -548 -397 -102 321 234 750 66 -266 -544 629 -173 -36 579 -1271 1002 290 -572 -509 541 276 178 12 113 -755 -330 518 387 -31 -889 348 1038 -1318 885 1137 -765 -1588 1119 -3546 3219 300 570 -146 57 -388 -243 456 -1589 -641 -714 -154 1702 1801 -299 -333 50 22 -1126 1101 -249 -237 265 -898 711 443 -2400 3028 -413 -380 -510 -213 932 363 -2773 489 652 886 370 -61 220 424 102 -394 -554 -2147 497 731 1573 -176 -471 178 -2324 1771 109 4 73 -194 500 -793 -1235 1088 -3120 2704 338 399 118 -212 72 -120 -277 386 170 -96 -1176 435 699 -350 769 -291 -2378 1236 339 271 -143 -755 -1158 959 -1071 1297 276 3 796 -953 699 114 213 -326 -129 326 -1217 -1889 2617 298 -430 -447 567 213 217 -31 -132 -838 -308 972 347 -1090 408 -859 -1198 -257 563 -155 1140 121 576 -514 -191 -343 285 132 904 826 423 -1096 -1049 1222 -1421 1920 -883 -721 -182 -13 637 441 -84 -129 -1862 676 455 1907 -69 -195 -619 567 -780 381 885 70 -818 390 -1319 115 837 591 -1573 1028 1 -1662 -640 1893 500 247 -1164 1149 164 -1629 440 648 996 289 -1575 -755 657 933 41 -2529 3156 -1796 872 -643 -90 623 960 466 -813 -451 793 387 549 465 -120 -1814 3458 -1944 261 -60 137 641 1491 -4990 1624 562 -456 596 435 -674 776 -194 -105 836 -336 21 3 397 -1688 416 333 -1056 547 989 1062 -328 -762 321 -2324 2198 725 -666 -5597 2187 2419 563 -1170 1207 -2013 2084 -732 168 429 252 -125 587 -524 -86 -545 -1011 1855 458 399 -1311 -3275 2272 2309 -40 255 655 -1165 -2770 1714 -813 1510 -36 -198 55 -548 -300 1112 483 -3211 3993 -479 -643 3 54 439 242 -238 -900 66 -1692 2449 -58 302 37 -724 -1879 1726 -614 -306 -82 1635 -567 -112 663 -107 -456 -1416 224 -156 696 433 743 -1162 -34 -1367 -168 1897 -279 974 182 53 -530 -938 799 -44 -1073 878 -721 -1411 2951 -444 626 581 -1312 43 -219 -88 508 111 364 -319 -86 -351 -1360 1549 260 261 -248 -1500 513 245 339 -99 563 258 -457 717 -3316 1981 476 369 390 -183 -1529 -798 1662 123 657 16 -363 -1210 -330 224 830 248 -909 1301 -1528 2251 -702 242 103 -66 205 705 -1381 -464 -2796 3518 129 447 228 -488 -471 102 195 -340 22 1140 -1666 -111 -2139 2933 -244 828 -191 -4455 1968 914 1299 -121 -288 -173 -470 -764 1659 -73 48 -2085 2666 -1266 234 408 227 -334 85 408 -3393 -4437 1074 4470 420 -139 -709 -31 152 427 -651 280 677 544 316 -583 -2175 1401 -50 253 -69 -960 830 135 -488 -2721 1485 -1633 147 2663 -1128 1301 63 345 -477 -542 1585 372 -877 -354 -367 574 -2354 1942 331 -182 213 79 -564 -1261 799 972 -290 -1139 -505 -1874 38 -867 93 -572 1673 981 -1754 455 933 

fore_arima_lt_2 = fore.arima.wge(daily_bike_data$Total_Users, phi = est$phi, theta = est$theta, d = 1, n.ahead = 60)
y.arma -184 548 213 38 6 -96 -551 -137 499 -58 -101 244 15 -173 -44 -204 -317 967 277 -384 -562 5 430 569 -1479 -75 736 -69 -2 405 -141 166 24 158 -703 618 89 -182 75 -67 208 -274 117 324 -98 300 360 452 -1292 177 -705 343 467 -110 -346 508 433 -956 405 283 -449 259 133 -1472 1267 261 -242 -1268 1354 155 285 -371 10 136 552 495 -122 -646 -394 626 -582 -256 345 286 -803 335 397 -889 149 542 25 997 -134 -1320 1013 333 -1670 984 440 453 -1314 128 1105 -141 -2331 2949 -315 -225 740 245 -2506 2353 155 -118 327 -528 186 537 717 -1961 1050 50 -1818 1800 175 106 -381 29 441 -621 682 -759 -696 1144 -595 165 -268 720 342 888 -1145 -386 218 486 -301 2 79 30 -690 -116 -8 994 344 30 -436 -358 285 -432 -486 671 380 -506 560 -129 289 -1413 1077 275 -375 -734 825 -328 283 201 211 103 -597 -60 577 290 -153 -243 -470 1394 -1378 -36 -37 -552 1296 -455 -795 172 84 742 454 385 -621 -844 83 -209 -548 -397 -102 321 234 750 66 -266 -544 629 -173 -36 579 -1271 1002 290 -572 -509 541 276 178 12 113 -755 -330 518 387 -31 -889 348 1038 -1318 885 1137 -765 -1588 1119 -3546 3219 300 570 -146 57 -388 -243 456 -1589 -641 -714 -154 1702 1801 -299 -333 50 22 -1126 1101 -249 -237 265 -898 711 443 -2400 3028 -413 -380 -510 -213 932 363 -2773 489 652 886 370 -61 220 424 102 -394 -554 -2147 497 731 1573 -176 -471 178 -2324 1771 109 4 73 -194 500 -793 -1235 1088 -3120 2704 338 399 118 -212 72 -120 -277 386 170 -96 -1176 435 699 -350 769 -291 -2378 1236 339 271 -143 -755 -1158 959 -1071 1297 276 3 796 -953 699 114 213 -326 -129 326 -1217 -1889 2617 298 -430 -447 567 213 217 -31 -132 -838 -308 972 347 -1090 408 -859 -1198 -257 563 -155 1140 121 576 -514 -191 -343 285 132 904 826 423 -1096 -1049 1222 -1421 1920 -883 -721 -182 -13 637 441 -84 -129 -1862 676 455 1907 -69 -195 -619 567 -780 381 885 70 -818 390 -1319 115 837 591 -1573 1028 1 -1662 -640 1893 500 247 -1164 1149 164 -1629 440 648 996 289 -1575 -755 657 933 41 -2529 3156 -1796 872 -643 -90 623 960 466 -813 -451 793 387 549 465 -120 -1814 3458 -1944 261 -60 137 641 1491 -4990 1624 562 -456 596 435 -674 776 -194 -105 836 -336 21 3 397 -1688 416 333 -1056 547 989 1062 -328 -762 321 -2324 2198 725 -666 -5597 2187 2419 563 -1170 1207 -2013 2084 -732 168 429 252 -125 587 -524 -86 -545 -1011 1855 458 399 -1311 -3275 2272 2309 -40 255 655 -1165 -2770 1714 -813 1510 -36 -198 55 -548 -300 1112 483 -3211 3993 -479 -643 3 54 439 242 -238 -900 66 -1692 2449 -58 302 37 -724 -1879 1726 -614 -306 -82 1635 -567 -112 663 -107 -456 -1416 224 -156 696 433 743 -1162 -34 -1367 -168 1897 -279 974 182 53 -530 -938 799 -44 -1073 878 -721 -1411 2951 -444 626 581 -1312 43 -219 -88 508 111 364 -319 -86 -351 -1360 1549 260 261 -248 -1500 513 245 339 -99 563 258 -457 717 -3316 1981 476 369 390 -183 -1529 -798 1662 123 657 16 -363 -1210 -330 224 830 248 -909 1301 -1528 2251 -702 242 103 -66 205 705 -1381 -464 -2796 3518 129 447 228 -488 -471 102 195 -340 22 1140 -1666 -111 -2139 2933 -244 828 -191 -4455 1968 914 1299 -121 -288 -173 -470 -764 1659 -73 48 -2085 2666 -1266 234 408 227 -334 85 408 -3393 -4437 1074 4470 420 -139 -709 -31 152 427 -651 280 677 544 316 -583 -2175 1401 -50 253 -69 -960 830 135 -488 -2721 1485 -1633 147 2663 -1128 1301 63 345 -477 -542 1585 372 -877 -354 -367 574 -2354 1942 331 -182 213 79 -564 -1261 799 972 -290 -1139 -505 -1874 38 -867 93 -572 1673 981 -1754 455 933 

t = 1:800
plot(t[670:731],daily_bike_data$Total_Users[670:731], type = 'l', main = "Short Term Forecast", xlim = c(670,795), xlab = "Time", ylab = "Total Users")
points(t[732:738],fore_arima_st_2$f, type = 'l', col = 'blue')
points(t[732:738],fore_arima_st_2$ll, type = 'l',lwd=2, lty = 2, col = 'red')
points(t[732:738],fore_arima_st_2$ul, type = 'l',lwd=2, lty = 2, col = 'red')

plot(t[670:731],daily_bike_data$Total_Users[670:731], type = 'l', main = "Long Term Forecast", xlim = c(670,795), xlab = "Time", ylab = "Total Users")
points(t[732:791],fore_arima_lt_2$f, type = 'l', col = 'blue')
points(t[732:791],fore_arima_lt_2$ll, type = 'l',lwd=2, lty = 2, col = 'red')
points(t[732:791],fore_arima_lt_2$ul, type = 'l',lwd=2, lty = 2, col = 'red')

  1. The models in factored form or at least separate the stationary and non- stationary factors with standard deviation or variance of the white noise.
  2. AIC
  3. ASE (short and long term forecasts)
  4. Rolling Window RMSE (short and long term forecasts)
  5. At least 10 superimposed spectral densities from 10 generated realizations like we did in Unit 11. Use these to help choose between the at least two candidate univariate models.
  6. Visualization of Forecasts for both the short- and long-term Horizons.
  7. Be sure and include confidence intervals when possible (I don’t have code for confidence intervals from MLP models at the moment… but that would be a good thing to work on!  )

Multivariate: At least one multivariate model (VAR or MLR with Correlated Errors)

  1. Include an ASE (rolling window is not yet available in multivariate models)
  1. Short Horizon (you pick the length.. could be one step ahead)
  2. Long Horizon (you pick … just must be longer than the short horizon.)
  1. Describe the explanatory variable(s) used in the model and why you felt they were significant / important.
  2. Visualization of Forecasts for both the short- and long-Horizons.
  3. Be sure and include confidence intervals if using VAR …
ccf(daily_bike_data$Total_Users,daily_bike_data$Humidity) # lag of 5 

daily_bike_data$Humidity_lag = dplyr::lag(daily_bike_data$Humidity,5)

fit = lm(Total_Users~Humidity_lag + Temperature + Temperature_Feels + Day_of_the_Week + Wind_Speed, data = daily_bike_data)

plotts.sample.wge(fit$residuals)$xbar

[1] 9.503264e-15
aic5.wge(fit$residuals, p =0:10, q = 0:5, type = 'bic')
---------WORKING... PLEASE WAIT... 


Five Smallest Values of  bic 
    p    q        bic
    2    1   13.57138
    1    3   13.57973
    3    1   13.58033
    2    2   13.58080
    4    1   13.58595
mlr = arima(daily_bike_data$Total_Users, order = c(6,0,1),xreg = cbind(daily_bike_data$Humidity_lag, daily_bike_data$Temperature, daily_bike_data$Temperature_Feels, daily_bike_data$Day_of_the_Week, daily_bike_data$Wind_Speed))
Warning in log(s2): NaNs produced
Warning in arima(daily_bike_data$Total_Users, order = c(6, 0, 1), xreg =
cbind(daily_bike_data$Humidity_lag, : possible convergence problem: optim gave
code = 1
plotts.wge(mlr$residuals) # looks random

acf(mlr$residuals[-(1:5)]) # only 1/20 acfs out of bounds 

ljung.wge(mlr$residuals) # greater than 0.05
Obs -0.01225566 -0.02716112 -0.03623396 -0.04795925 -0.05746748 -0.0794555 -0.02938715 -0.001984348 -0.01131826 0.05093032 -0.003201458 -0.002334059 -0.02369733 0.03374261 0.06930404 -0.01419825 0.02123568 -0.01049807 -0.007215088 0.01094336 0.02458074 -0.01208151 -0.02100042 0.001923981 
$test
[1] "Ljung-Box test"

$K
[1] 24

$chi.square
[1] 19.56278

$df
[1] 24

$pval
[1] 0.7213486
ljung.wge(mlr$residuals, K =48)
Obs -0.01225566 -0.02716112 -0.03623396 -0.04795925 -0.05746748 -0.0794555 -0.02938715 -0.001984348 -0.01131826 0.05093032 -0.003201458 -0.002334059 -0.02369733 0.03374261 0.06930404 -0.01419825 0.02123568 -0.01049807 -0.007215088 0.01094336 0.02458074 -0.01208151 -0.02100042 0.001923981 -0.02056452 0.05950971 0.0430777 0.04935224 0.1186554 -0.04045242 -0.01594315 -0.03858605 0.02273316 0.02683504 -0.02047584 0.01395783 0.09418566 -0.04322951 0.04973698 -0.02609069 0.05379182 -0.01609637 -0.03581977 -0.002512819 -0.03853339 0.04890081 0.02738517 -0.007969906 
$test
[1] "Ljung-Box test"

$K
[1] 48

$chi.square
[1] 58.45805

$df
[1] 48

$pval
[1] 0.1433528
mlr_st_pred = predict(mlr, newxreg = data.frame(daily_bike_data$Humidity_lag[725:731], daily_bike_data$Temperature[725:731], daily_bike_data$Temperature_Feels[725:731], daily_bike_data$Day_of_the_Week[725:731], daily_bike_data$Wind_Speed[725:731]), n.ahead = 7, lastn= TRUE)

plot(seq(1,731,1),daily_bike_data$Total_Users,type = 'l',xlim = c(720,731))
points(seq(725,731,1),mlr_st_pred$pred,type = 'b', pch = 15)

mlr_st_ASE= mean((daily_bike_data$Total_Users[725:731] - mlr_st_pred$pred)^2)

mlr_st_ASE
[1] 1287463
mlr_lt_pred = predict(mlr, newxreg = data.frame(daily_bike_data$Humidity_lag[672:731], daily_bike_data$Temperature[672:731], daily_bike_data$Temperature_Feels[672:731], daily_bike_data$Day_of_the_Week[672:731], daily_bike_data$Wind_Speed[672:731]), n.ahead = 60, lastn= TRUE)

plot(seq(1,731,1),daily_bike_data$Total_Users,type = 'l',xlim = c(670,731))
points(seq(672,731,1),mlr_lt_pred$pred,type = 'b', pch = 15)

mlr_lt_ASE= mean((daily_bike_data$Total_Users[672:731] - mlr_lt_pred$pred)^2)

mlr_lt_ASE
[1] 3731991
mlr_st_pred2 = predict(mlr, newxreg = data.frame(daily_bike_data$Humidity_lag[725:731], daily_bike_data$Temperature[725:731], daily_bike_data$Temperature_Feels[725:731], daily_bike_data$Day_of_the_Week[725:731], daily_bike_data$Wind_Speed[725:731]), n.ahead = 7, lastn= FALSE)
plot(seq(1,731,1),daily_bike_data$Total_Users,type = 'l',xlim = c(720,738), main = "Short Term MLR Forecast")
points(seq(732,738,1),mlr_st_pred2$pred,type = 'b', pch = 15,col = 'blue')

mlr_lt_pred2 = predict(mlr, newxreg = data.frame(daily_bike_data$Humidity_lag[672:731], daily_bike_data$Temperature[672:731], daily_bike_data$Temperature_Feels[672:731], daily_bike_data$Day_of_the_Week[672:731], daily_bike_data$Wind_Speed[672:731]), n.ahead = 60, lastn= FALSE)
plot(seq(1,731,1),daily_bike_data$Total_Users,type = 'l',xlim = c(720,785), main = "Long Term MLR Forecast")
points(seq(732,791,1),mlr_lt_pred2$pred,type = 'b', pch = 15,col = 'blue')

VARselect(daily_bike_data[,c(5,8:11,14)])
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
     7      7      7      7 

$criteria
                  1            2            3            4            5
AIC(n) 2.762319e+01 2.733958e+01 2.716829e+01 2.689074e+01 2.632813e+01
HQ(n)  2.772619e+01 2.753087e+01 2.744788e+01 2.725862e+01 2.678429e+01
SC(n)  2.789002e+01 2.783512e+01 2.789256e+01 2.784372e+01 2.750982e+01
FPE(n) 9.922004e+11 7.472005e+11 6.296103e+11 4.770590e+11 2.718306e+11
                   6             7             8             9            10
AIC(n) -3.750199e+01 -4.271176e+01 -4.266486e+01 -4.261632e+01 -4.256334e+01
HQ(n)  -3.695753e+01 -4.207901e+01 -4.194382e+01 -4.180699e+01 -4.166572e+01
SC(n)  -3.609158e+01 -4.107264e+01 -4.079703e+01 -4.051977e+01 -4.023808e+01
FPE(n)  5.168071e-17  2.824148e-19  2.960966e-19  3.109863e-19  3.281202e-19
fit2 = VAR(daily_bike_data[,c(5,8:11,14)], p = 10, type = 'both')


plotts.wge(fit2$varresult$Total_Users$residuals)

acf(fit2$varresult$Total_Users$residuals)

ljung.wge(fit2$varresult$Total_Users$residuals, p = 10)
Obs -0.0005959342 -0.004132474 0.005197177 -0.0098346 -0.005829737 -0.01447498 -0.005342387 -0.02915348 -0.03548842 -0.02164495 -0.0445405 -0.03091929 -0.05910042 0.0379835 0.06817437 -0.03618838 0.01217303 -0.0167584 -0.0230686 -0.0006933888 0.04627761 0.009535346 -0.003221014 0.01358627 
$test
[1] "Ljung-Box test"

$K
[1] 24

$chi.square
[1] 14.90476

$df
[1] 14

$pval
[1] 0.3846962
ljung.wge(fit2$varresult$Total_Users$residuals, p = 10, K = 48)
Obs -0.0005959342 -0.004132474 0.005197177 -0.0098346 -0.005829737 -0.01447498 -0.005342387 -0.02915348 -0.03548842 -0.02164495 -0.0445405 -0.03091929 -0.05910042 0.0379835 0.06817437 -0.03618838 0.01217303 -0.0167584 -0.0230686 -0.0006933888 0.04627761 0.009535346 -0.003221014 0.01358627 -0.0301048 0.04538141 0.01339568 0.03353197 0.1077264 -0.04467969 0.01035195 -0.05387493 0.009563282 0.01698439 0.002288517 -0.02275094 0.05485335 -0.06199872 0.03351602 -0.03592237 0.03521859 -0.003384203 -0.01519705 0.01945388 -0.03743808 0.01342283 0.02967684 0.005916222 
$test
[1] "Ljung-Box test"

$K
[1] 48

$chi.square
[1] 41.74985

$df
[1] 38

$pval
[1] 0.3111088
var_st_pred = predict(fit2, n.ahead = 7, lastn = TRUE)
Warning in summary.lm(x): essentially perfect fit: summary may be unreliable
t = 1:731
plot(seq(1,731,1),daily_bike_data$Total_Users,type = 'l',xlim = c(670,731))
points(t[725:731], var_st_pred$fcst$Total_Users[,1], type="l", lwd=2, lty = 1)

var_st_ase = mean((daily_bike_data$Total_Users[725:731]-var_st_pred$fcst$Total_Users[,1])^2)
var_st_ase
[1] NA
var_lt_pred = predict(fit2, n.ahead = 60, lastn = TRUE)
Warning in summary.lm(x): essentially perfect fit: summary may be unreliable
t = 1:731
plot(seq(1,731,1),daily_bike_data$Total_Users,type = 'l',xlim = c(672,731))
points(t[672:731], var_lt_pred$fcst$Total_Users[,1], type="l", lwd=2, lty = 1)

var_lt_ase = mean((daily_bike_data$Total_Users[672:731]-var_lt_pred$fcst$Total_Users[,1])^2)
var_lt_ase
[1] NA
var_st_pred2 = predict(fit2, n.ahead = 7, lastn = FALSE)
Warning in summary.lm(x): essentially perfect fit: summary may be unreliable
t = 1:800
plot(seq(1,731,1),daily_bike_data$Total_Users,type = 'l',xlim = c(670,738), main = "Short Term VAR Forecast")
points(t[732:738], var_st_pred2$fcst$Total_Users[,1], type="l", lwd=2, lty = 1, col = 'blue')

var_lt_pred2 = predict(fit2, n.ahead = 60, lastn = FALSE)
Warning in summary.lm(x): essentially perfect fit: summary may be unreliable
t = 1:800
plot(seq(1,731,1),daily_bike_data$Total_Users,type = 'l',xlim = c(670,791), main = "Long Term VAR Forecast")
points(t[732:791], var_lt_pred2$fcst$Total_Users[,1], type="l", lwd=2, lty = 1, col = 'blue')

MLP Model

  1. ASE (short and long term forecasts)
  2. Rolling Window RMSE (short and long term forecasts (only if univariate)
  3. Visualization of Forecasts for both the short- and long-term Horizons.
  4. Confidence / Prediction intervals are not required (I don’t have code for confidence / prediction intervals (bootstrap intervals) for MLP models at the moment… but that would be a good thing to work on!  )
# Convert tibble to a data frame with ts() objects
bike_DF <- daily_bike_data %>%
  dplyr::select(-c(Date, Hour, Holiday, Temperature_Feels, Casual_Users, Registered_Users)) %>%  # Remove unnecessary columns
  mutate(across(everything(), ~ zoo(.x, order.by = daily_bike_data$Date))) %>%  # Convert each column to a zoo object
  mutate(across(everything(), ~ as.ts(.x))) %>%  # Convert zoo objects to ts objects
  as.data.frame()  # Convert the tibble to a data frame

bikeShortTrain = bike_DF[1:724,]
bikeShortTest = bike_DF[725:731,]
bikeLongTrain = bike_DF[1:671,]
bikeLongTest = bike_DF[672:731,]
seed = 137
# Forecast Short-term predictor variables using MLP
set.seed(seed)
fit.mlp.short.Season = mlp(ts(bikeShortTrain[,"Season"]),reps = 10, difforder = 0, comb = "mean", det.type = "bin")
fore.mlp.short.Season = forecast(fit.mlp.short.Season, h = 7)
plot(fore.mlp.short.Season)

fit.mlp.short.Day_of_the_Week = mlp(ts(bikeShortTrain[,"Day_of_the_Week"]),reps = 10, difforder = 0, comb = "mean", det.type = "bin")
fore.mlp.short.Day_of_the_Week = forecast(fit.mlp.short.Day_of_the_Week, h = 7)
fore.mlp.short.Day_of_the_Week$mean = round(fore.mlp.short.Day_of_the_Week$mean) # Round to nearest whole number
plot(fore.mlp.short.Day_of_the_Week)

fit.mlp.short.Working_Day = mlp(ts(bikeShortTrain[,"Working_Day"]),reps = 10, difforder = 0, comb = "mean", det.type = "bin")
fore.mlp.short.Working_Day = forecast(fit.mlp.short.Working_Day, h = 7)
fore.mlp.short.Working_Day$mean = round(fore.mlp.short.Working_Day$mean) # Round to nearest whole number
plot(fore.mlp.short.Working_Day)

# fit.mlp.short.Weather_Type = mlp(ts(bikeShortTrain[,"Weather_Type"]),reps = 5, difforder = 0, comb = "mean", det.type = "auto")
# fore.mlp.short.Weather_Type = forecast(fit.mlp.short.Weather_Type, h = 7)
# plot(fore.mlp.short.Weather_Type)

fit.mlp.short.Temperature = mlp(ts(bikeShortTrain[,"Temperature"]),reps = 10, difforder = 0, comb = "mean", det.type = "bin")
fore.mlp.short.Temperature = forecast(fit.mlp.short.Temperature, h = 7)
plot(fore.mlp.short.Temperature)

fit.mlp.short.Humidity = mlp(ts(bikeShortTrain[,"Humidity"]),reps = 10, difforder = 0, comb = "mean", det.type = "bin")
fore.mlp.short.Humidity = forecast(fit.mlp.short.Humidity, h = 7)
plot(fore.mlp.short.Humidity)

fit.mlp.short.Wind_Speed = mlp(ts(bikeShortTrain[,"Wind_Speed"]),reps = 10, difforder = 0, comb = "mean", det.type = "bin")
fore.mlp.short.Wind_Speed = forecast(fit.mlp.short.Wind_Speed, h = 7)
plot(fore.mlp.short.Wind_Speed)

# Forecast Long-term predictor variables using MLP
set.seed(seed)
fit.mlp.long.Season = mlp(ts(bikeLongTrain[,"Season"]),reps = 10, difforder = 0, comb = "median", det.type = "bin")
fore.mlp.long.Season = forecast(fit.mlp.long.Season, h = 60)
plot(fore.mlp.long.Season)

fit.mlp.long.Day_of_the_Week = mlp(ts(bikeLongTrain[,"Day_of_the_Week"]),reps = 10, difforder = 0, comb = "median", det.type = "bin")
fore.mlp.long.Day_of_the_Week = forecast(fit.mlp.long.Day_of_the_Week, h = 60)
fore.mlp.long.Day_of_the_Week$mean = round(fore.mlp.long.Day_of_the_Week$mean) # Round to nearest whole number
plot(fore.mlp.long.Day_of_the_Week)

fit.mlp.long.Working_Day = mlp(ts(bikeLongTrain[,"Working_Day"]),reps = 10, difforder = 0, comb = "median", det.type = "bin")
fore.mlp.long.Working_Day = forecast(fit.mlp.long.Working_Day, h = 60)
fore.mlp.long.Working_Day$mean = round(fore.mlp.long.Working_Day$mean) # Round to nearest whole number
plot(fore.mlp.long.Working_Day)

# fit.mlp.long.Weather_Type = mlp(ts(bikeLongTrain[,"Weather_Type"]),reps = 10, difforder = 0, comb = "mean", det.type = "auto")
# fore.mlp.long.Weather_Type = forecast(fit.mlp.long.Weather_Type, h = 60)
# plot(fore.mlp.long.Weather_Type)

fit.mlp.long.Temperature = mlp(ts(bikeLongTrain[,"Temperature"]),reps = 10, difforder = 0, comb = "mean", det.type = "bin")
fore.mlp.long.Temperature = forecast(fit.mlp.long.Temperature, h = 60)
plot(fore.mlp.long.Temperature)

fit.mlp.long.Humidity = mlp(ts(bikeLongTrain[,"Humidity"]),reps = 10, difforder = 0, comb = "mean", det.type = "bin")
fore.mlp.long.Humidity = forecast(fit.mlp.long.Humidity, h = 60)
plot(fore.mlp.long.Humidity)

fit.mlp.long.Wind_Speed = mlp(ts(bikeLongTrain[,"Wind_Speed"]),reps = 10, difforder = 0, comb = "mean", det.type = "bin")
fore.mlp.long.Wind_Speed = forecast(fit.mlp.long.Wind_Speed, h = 60)
plot(fore.mlp.long.Wind_Speed)

# Pack into data frames
BDF_fore_short = data.frame(Day = ts(seq(1,731,1)),
                            Season = ts(c(bikeShortTrain$Season,fore.mlp.short.Season$mean)),
                            Day_of_the_Week = ts(c(bikeShortTrain$Day_of_the_Week,fore.mlp.short.Day_of_the_Week$mean)),
                            Working_Day = ts(c(bikeShortTrain$Working_Day,fore.mlp.short.Working_Day$mean)),
                            # Weather_Type = ts(c(bikeShortTrain$Weather_Type,fore.mlp.short.Weather_Type$mean)),
                            Temperature = ts(c(bikeShortTrain$Temperature,fore.mlp.short.Temperature$mean)),
                            Humidity = ts(c(bikeShortTrain$Humidity,fore.mlp.short.Humidity$mean)),
                            Wind_Speed = ts(c(bikeShortTrain$Wind_Speed,fore.mlp.short.Wind_Speed$mean)))
BDF_fore_long = data.frame(Day = ts(seq(1,731,1)),
                            Season = ts(c(bikeLongTrain$Season,fore.mlp.long.Season$mean)),
                            Day_of_the_Week = ts(c(bikeLongTrain$Day_of_the_Week,fore.mlp.long.Day_of_the_Week$mean)),
                            Working_Day = ts(c(bikeLongTrain$Working_Day,fore.mlp.long.Working_Day$mean)),
                            # Weather_Type = ts(c(bikeLongTrain$Weather_Type,fore.mlp.long.Weather_Type$mean)),
                            Temperature = ts(c(bikeLongTrain$Temperature,fore.mlp.long.Temperature$mean)),
                            Humidity = ts(c(bikeLongTrain$Humidity,fore.mlp.long.Humidity$mean)),
                            Wind_Speed = ts(c(bikeLongTrain$Wind_Speed,fore.mlp.long.Wind_Speed$mean)))

BDF_xreg_short = data.frame(Day = ts(seq(1,724,1)),
                            Season = ts(c(bikeShortTrain$Season)),
                            Day_of_the_Week = ts(c(bikeShortTrain$Day_of_the_Week)),
                            Working_Day = ts(c(bikeShortTrain$Working_Day)),
                            # Weather_Type = ts(c(bikeShortTrain$Weather_Type)),
                            Temperature = ts(c(bikeShortTrain$Temperature)),
                            Humidity = ts(c(bikeShortTrain$Humidity)),
                            Wind_Speed = ts(c(bikeShortTrain$Wind_Speed)))
BDF_xreg_long = data.frame(Day = ts(seq(1,671,1)),
                            Season = ts(c(bikeLongTrain$Season)),
                            Day_of_the_Week = ts(c(bikeLongTrain$Day_of_the_Week)),
                            Working_Day = ts(c(bikeLongTrain$Working_Day)),
                            # Weather_Type = ts(c(bikeLongTrain$Weather_Type)),
                            Temperature = ts(c(bikeLongTrain$Temperature)),
                            Humidity = ts(c(bikeLongTrain$Humidity)),
                            Wind_Speed = ts(c(bikeLongTrain$Wind_Speed)))

# Fit MLP model for Total_Users
set.seed(seed)
fit.mlp.st = mlp(ts(bikeShortTrain$Total_Users),reps = 10, comb = "median",xreg = BDF_xreg_short, allow.det.season = TRUE, det.type = "bin")
fit.mlp.lt = mlp(ts(bikeLongTrain$Total_Users),reps = 10, comb = "median",xreg = BDF_xreg_long, allow.det.season = TRUE, det.type = "bin")

# Forecast and evaluate ASE for short and long
fore.mlp.st = forecast(fit.mlp.st, h = 7, xreg = BDF_fore_short)
fore.mlp.lt = forecast(fit.mlp.lt, h = 60, xreg = BDF_fore_long)
ASE.mlp.st = mean((bikeShortTest$Total_Users - fore.mlp.st$mean)^2)
ASE.mlp.lt = mean((bikeLongTest$Total_Users - fore.mlp.lt$mean)^2)
print(paste("MLP ASE Short Term: ", round(ASE.mlp.st, 2)))
[1] "MLP ASE Short Term:  4177344.82"
print(paste("MLP ASE Long Term: ", round(ASE.mlp.lt, 2)))
[1] "MLP ASE Long Term:  10000456.14"
plot(fore.mlp.st)

plot(fore.mlp.lt)

# Plot the forecasts
t = 1:731
plot(t[720:731],bike_DF$Total_Users[720:731], type = 'l', xlab = "Time", ylab = "Total Users", main = "MLP Short Term Forecast")
lines(t[725:731], fore.mlp.st$mean, type="l", lwd=2, lty = 2, col = 'blue')

plot(t[600:731],bike_DF$Total_Users[600:731], type = 'l', xlab = "Time", ylab = "Total Users", main = "MLP Long Term Forecast")
lines(t[672:731], fore.mlp.lt$mean, type="l", lwd=2, lty = 2, col = 'blue')

Ensemble Model

# Code for Figure 11.37
# TODO: Adopt this for ours
#Ensemble

#VAR p = 7 non seasonal
CMortVAR7 = VAR(cardiacTrain, type = "both", p = 7) #p = 2 from SBC
preds7=predict(CMortVAR7,n.ahead=156)
RMSEVAR7 = sqrt(mean((cardiacTest[,"cmort"] - preds7$fcst$cmort[,1])^2))
RMSEVAR7 # 6.664

#VAR p = 2 seasonal
CMortVAR2S = VAR(cardiacTrain, season = 52, type = "both", p = 2) #p = 2 from SBC
preds2S=predict(CMortVAR2S,n.ahead=156)
ensemble = (preds2S$fcst$cmort[,1] + fore.mlp.cmort$mean)/2

#Plot
plot(seq(1,508,1), cardiac[,"cmort"], type = "l",xlim = c(350,508), ylim = c(70,110), xlab = "Time", ylab = "Cardiac Mortality", main = "52 Week Cardiac Mortality Forecast From A VAR/MLP Ensemble")
lines(seq(353,508,1), ensemble, type = "l", lwd = 4, col = "green")
lines(seq(353,508,1),preds2S$fcst$cmort[,1] , type = "l", lwd = 2, lty = 2, col = "red")
lines(seq(353,508,1),fore.mlp.cmort$mean , type = "l", lwd = 2, lty = 4, col = "blue")
lines(seq(353,508,1),preds7$fcst$cmort[,1] , type = "l", lwd = 2, lty = 2, col = "purple")
RMSEVAR7 = sqrt(mean((cardiacTest[,"cmort"] - preds7$fcst$cmort[,1])^2))
RMSEVAR7
RMSEVAR2S = sqrt(mean((cardiacTest[,"cmort"] - preds2S$fcst$cmort[,1])^2))
RMSEVAR2S
RMSE = sqrt(mean((cardiacTest[,"cmort"] - fore.mlp.cmort$mean)^2))
RMSE
RMSEENSEMBLE = sqrt(mean((cardiacTest[,"cmort"] - ensemble)^2))
RMSEENSEMBLE # 5.64, 5,81

Model Comparison and Final Forecasts

  1. Provide a table comparing all models on at least ASE and rwRMSE (if available).
  2. Include at least one ensemble model in addition to the models above.
  3. Make a case as to why you feel one of your candidate models is the most useful.
  4. Provide you final short and long term forecasts with that model.
# Print a table with the ASE from all models:

Conclusion